#!/usr/bin/perl

# UPDATED SO THAT WE GET PRECISE START END OUTPUT

use strict;
use warnings;

sub trim {
my $x = shift;
$x =~ s/^\s*//g;
$x =~ s/\s*$//g;
return $x;
}
# PHASE 1 - detect valid complexes

my $threshold = 2.0;

open TABLE, "helix_table.txt";
#open TABLE, "ht2.txt";
my $lastpdb = "";
my $pdb = "";
my %helices;

open OUT, ">ccids_4.0_verbose_10A.txt";

use Class::Struct;
struct(helix  => [ddg => '$', start => '$', end => '$']);
struct(point  => [rn => '$', x => '$', y => '$', z => '$', ddg => '$']);

use Graph::Directed;
use Math::Trig;
# find PDBs in which 

foreach my $line (<TABLE>) {
	chomp $line;
	if ($line =~ /^$/) {last;}
	my @objs = split ('\t',$line);	
	# figure out all the chains this pdb code corresponds to
	$pdb = $objs[0];
	#print $objs[7];
	#print $objs[8];
	#last;
	if ($pdb ne $lastpdb) {
		# end prior pdb's network
		analyze (%helices);

		undef %helices;
	} 

	# push total ddg to array of helices between chains in question
	#push ($helices{$objs[8]}{(substr $objs[7], 2, 1)}, $objs[11]); 
	#if (defined ($helices{$objs[8]}{(substr $objs[7], 2, 1)}) {
		my $a = new helix;
		$a -> ddg  ($objs[11]);
		$a -> start ($objs[23]); 
		$a -> end   ($objs[24]); 
		push (@{$helices{$objs[8]}{(substr $objs[7], 2, 1)}}, $a);
	#} else {
	#}
	#-> ddg   ($objs[11]); 
	#$helices{$objs[8]}{(substr $objs[7], 2, 1)}-> start ($objs[23]); 
	#$helices{$objs[8]}{(substr $objs[7], 2, 1)}-> end   ($objs[24]); 
	$lastpdb = $pdb;
}


close OUT;

sub analyze {
	my (%helices) = @_;
	my $heligraph = Graph::Directed->new;
	print "$lastpdb\n"; 
	my $interesting = 0;
	foreach my $ch1 (keys %helices) {
		#print "$ch1 has a helix interacting with ";
		foreach my $ch2 (keys %{$helices{$ch1}}) {
			#print "$ch2 for " . $helices{$ch1}{$ch2} . ", ";
			if (scalar(@{$helices{$ch1}{$ch2}}) > 1) {
				my @starts;
				my @ends;
				foreach my $h (@{$helices{$ch1}{$ch2}}) {
					if ($h->ddg < $threshold) { next; }
					
					open PDB, "../geometry/pdb_partner_files/$lastpdb\_$ch1\_$ch2.pdb" or return;
					# record start and end of each and calc angle in between vectors
					my $s = new point; 
					my $e = new point;
					my $efound = 0;
					my $sfound = 0;
					while (<PDB>) {
						if (/^ATOM/) {
							if ((substr $_, 21, 1) eq $ch1 ) {
								if (trim(substr($_,22,4)) == $h->start) {	
									$s->x (trim(substr($_,30,8 ))); $s->y (trim(substr($_, 38,8)) ); $s->z (trim(substr($_,46,8)) );
									$s->rn ($h->start);
									$s->ddg ($h->ddg);
									$sfound = 1;
								} elsif (trim(substr($_,22,4)) == $h->end) {
									$e->x (trim(substr($_,30,8 ))); $e->y (trim(substr($_, 38,8)) ); $e->z (trim(substr($_,46,8)) );
									$e->rn ($h->end);
									$e->ddg ($h->ddg);
									$efound = 1;
								}
							}
						}
					}
					if ($efound == 0) {print "weird, didn't find end to $ch1 " . $h->start ." " . $h->end . " \n";}
					if ($sfound == 0) {print "weird, didn't find start to $ch1 " . $h->start ." " . $h->end . " \n";}
					push (@starts, $s);
					push (@ends, $e);
				}
				for (my $i = 0; $i < scalar(@starts); $i++) {
					my $vec1 = new point;
					$vec1->x ( ($ends[$i]->x - $starts[$i]->x));
					$vec1->y ( ($ends[$i]->y - $starts[$i]->y));
					$vec1->z ( ($ends[$i]->z - $starts[$i]->z));
					for (my $j = 1; $j < scalar(@starts); $j++) {
						my $vec2 = new point;
						$vec2->x ( ($ends[$j]->x - $starts[$j]->x));
						$vec2->y ( ($ends[$j]->y - $starts[$j]->y));
						$vec2->z ( ($ends[$j]->z - $starts[$j]->z));
						my $ang = rad2deg( acos( ( $vec1->x * $vec2->x + $vec1->y * $vec2->y + $vec1->z * $vec2->z ) / ( (sqrt($vec1->x**2 + $vec1->y**2 + $vec1->z**2) * (sqrt($vec2->x**2 + $vec2->y**2 + $vec2->z**2) ) ) ) ));
						#if ($ang < 0) {$ang = -$ang;}
						#if (!(180 - $ang < 20 || $ang < 5)) { next;}

						my $dist = new point;
						$dist->x ( ($starts[$j]->x - $starts[$i]->x));
						$dist->y ( ($starts[$j]->y - $starts[$i]->y));
						$dist->z ( ($starts[$j]->z - $starts[$i]->z));
						# mostly just don't want this parallel to helix axis
						# indication that it's two fragmnets of one helix
						if ($dist->x == 0 && $dist->y == 0 && $dist->z == 0) {next;}
						my $perang = rad2deg(acos( ( $vec1->x * $dist->x + $vec1->y * $dist->y + $vec1->z * $dist->z ) / ( (sqrt($vec1->x**2 + $vec1->y**2 + $vec1->z**2) * (sqrt($dist->x**2 + $dist->y**2 + $dist->z**2) ) ) ) ));
						#if ($perang > -10 && $perang < 10 || $perang < -170 && $perang > -190 || $perang > 170 && $perang < 190 || $perang > 350 || $perang < -350) {next;}
						
						print OUT "$lastpdb mult $ch1 " . $starts[$i]->rn . " " . $ends[$i]->rn . " ddg " . $starts[$i]->ddg . "  $ch1 " . $starts[$j]->rn . " " . $ends[$j]->rn . " ddg " . $starts[$j]->ddg . " ang $ang perang $perang vs $ch2\n";
					}
				}
			}
			foreach my $ch3 (keys %helices) {
				if (defined $helices{$ch3}{$ch2} ){
					if ( $ch1 eq $ch3 or $ch1 eq $ch2 or $ch2 eq $ch3) {next;}
					my @h12starts;
					my @h12ends;	
					foreach my $h12 (@{$helices{$ch1}{$ch2}}) {
						if ($h12->ddg < $threshold) {
							next;
						}	
					
						open PDB, "../geometry/pdb_partner_files/$lastpdb\_$ch1\_$ch2.pdb" or return;
						# record start and end of each and calc angle in between vectors
						my $h12s = new point; 
						my $h12e = new point;
						my $h12efound = 0;
						my $h12sfound = 0;
						while (<PDB>) {
							if (/^ATOM/) {
								if ((substr $_, 21, 1) eq $ch1 ) {
									if (trim(substr($_,22,4)) == $h12->start) {	
										$h12s->x (trim(substr($_,30,8 ))); $h12s->y (trim(substr($_, 38,8)) ); $h12s->z (trim(substr($_,46,8)) );
										$h12s->rn ($h12->start);
										$h12s->ddg ($h12->ddg);
										$h12sfound = 1;
									} elsif (trim(substr($_,22,4)) == $h12->end) {
										$h12e->x (trim(substr($_,30,8 ))); $h12e->y (trim(substr($_, 38,8)) ); $h12e->z (trim(substr($_,46,8)) );
										$h12e->rn ($h12->end);
										$h12e->ddg ($h12->ddg);
										$h12efound = 1;
									}
								}
							}
						}
						if ($h12efound == 0) {print "weird, didn't find end to $ch1 " . $h12->start ." " . $h12->end . " \n";}
						if ($h12sfound == 0) {print "weird, didn't find start to $ch1 " . $h12->start ." " . $h12->end . " \n";}
						push (@h12starts, $h12s);
						push (@h12ends, $h12e);
						my @h32starts;
						my @h32ends;
						foreach my $h32 (@{$helices{$ch3}{$ch2}}) {
							if ($h32->ddg < $threshold) {
								next;
							}	
					
							open PDB, "../geometry/pdb_partner_files/$lastpdb\_$ch3\_$ch2.pdb" or return;
							# record start and end of each and calc angle in between vectors
							my $h32s = new point; 
							my $h32e = new point;
							my $h32efound = 0;
							my $h32sfound = 0;
							while (<PDB>) {
								if (/^ATOM/) {
									if ((substr $_, 21, 1) eq $ch3 ) {
										if (trim(substr($_,22,4)) == $h32->start) {	
											$h32s->x (trim(substr($_,30,8 ))); $h32s->y (trim(substr($_, 38,8)) ); $h32s->z (trim(substr($_,46,8)) );
											$h32s->rn ($h32->start);
											$h32s->ddg ($h32->ddg);
											$h32sfound = 1;
										} elsif (trim(substr($_,22,4)) == $h32->end) {
											$h32e->x (trim(substr($_,30,8 ))); $h32e->y (trim(substr($_, 38,8)) ); $h32e->z (trim(substr($_,46,8)) );
											$h32e->rn ($h32->end);
											$h32e->ddg ($h32->ddg);
											$h32efound = 1;
										}
									}
								}
							}
							if ($h32efound == 0) {print "weird, didn't find end to $ch3 " . $h32->start ." " . $h32->end . " \n";}
							if ($h32sfound == 0) {print "weird, didn't find start to $ch3 " . $h32->start ." " . $h32->end . " \n";}
							
							my $distsq = ($h32s->x - $h12s->x)**2 + ($h32s->y - $h12s->y)**2+ ($h32s->z - $h32s->z)**2;
							if ($distsq < 100 ) {
								push (@h32starts, $h32s);
								push (@h32ends, $h32e);
							}
							for (my $i = 0; $i < scalar(@h12starts); $i++) {
								my $vec1 = new point;
								$vec1->x ( ($h12ends[$i]->x - $h12starts[$i]->x));
								$vec1->y ( ($h12ends[$i]->y - $h12starts[$i]->y));
								$vec1->z ( ($h12ends[$i]->z - $h12starts[$i]->z));
								for (my $j = 0; $j < scalar(@h32starts); $j++) {
									my $vec2 = new point;
									$vec2->x ( ($h32ends[$j]->x - $h32starts[$j]->x));
									$vec2->y ( ($h32ends[$j]->y - $h32starts[$j]->y));
									$vec2->z ( ($h32ends[$j]->z - $h32starts[$j]->z));
									my $ang = rad2deg(acos( ( $vec1->x * $vec2->x + $vec1->y * $vec2->y + $vec1->z * $vec2->z ) / ( (sqrt($vec1->x**2 + $vec1->y**2 + $vec1->z**2) * (sqrt($vec2->x**2 + $vec2->y**2 + $vec2->z**2) ) ) ) ) );
									if ($ang < 0) {$ang = -$ang;}
									#if (180 - $ang < 5 || $ang < 5) {
										print OUT "$lastpdb coop $ch1 " . $h12starts[$i]->rn . " " . $h12ends[$i]->rn . " ddg " . $h12ends[$i]->ddg . " $ch3 " . $h32starts[$j]->rn . " " . $h32ends[$j]->rn . " ddg " . $h32ends[$j]->ddg . " ang $ang vs $ch2\n";
										print "$lastpdb coop $ch1 " . $h12starts[$i]->rn . " " . $h12ends[$i]->rn . " $ch3 " . $h32starts[$j]->rn . " " . $h32ends[$j]->rn . " vs $ch2\n";
									#}
								}
							}
						}
					}
				}
			}
		}
		#print "\n";
	}
}
